DISCONTINUOUS GALERKIN METHODS FOR THE 
j^BIHARMONIC EQUATION FROM A DISCRETE 
VARIATIONAL PERSPECTIVE 



TRISTAN PRYER 



Abstract. We study discontinuous Galerkin approximations of the p-biharmonic 
equation from a variational perspective. We propose a discrete variational for- 
mulation of the problem based on a appropriate definition of a finite element 
Hessian and study convergence of the method (writhout rates) using a weak 
lower semicontinuity argument. 

We present numerical experiments aimed at testing the robustness of the 
method. We also note a superconvergence effect for some values of p. 

1. Introduction, problem setup and notation 

The p-biharmonic equation is a fourth order elhptic boundary value problem, 
related to, in fact a nonlinear generalisation of, the biharmonic problem. Such 
problems typically arise from areas of elasticity. It is a fourth order analog to 
its second order sibling, the p-Laplacian, and as such is useful as a prototypical 
nonlinear fourth order problem. 

The efficient numerical simulation of general fourth order problems has attracted 
recent interest. A conforming approach to this class of problem would require 
the use of finite elements, the Argyris element for example |Cia781 §6]. From 
a practical point of view the approach presents difficulties, in that the finite 
elements are difficult to design and complicated to implement, especially when 
working in three spatial dimensions. 

Discontinuous Galerkin (dG) methods form a class of nonconforming finite ele- 
ment method. They are extremely popular due to their successful application to an 
ever expanding range of problems. A very accessible unification of these methods 
together with a detailed historical overview is presented in .ABCM02j . 

If p = 2 we have the special case that the (2-) biharmonic problem is linear. 
It has been well studied in the context of dG methods, for example, the papers 
|LS03[ [GH09j study the use of h~p dG finite elements (where p here means the 
local polynomial degree) applied to the (2-) biharmonic problem. To the authors 
knowledge there is currently no finite element method posed for the general p- 
biharmonic problem. 

In this work we use discrete variational techniques to build a discontinuous 
Galerkin (dG) numerical scheme for the p-biharmonic operator. We are interested 
in such a methodology due to the applications to discrete symmetries, in particular, 
discrete versions of Noether's Theorem |Noe71l rMP12aj . In a sibling publication 
|MP12bj we prove discrete versions of Noether's Theorem for general dG methods 
for second order variational minimisation problems which include the p-biharmonic 
problem as a prototypical example. 
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A key constituent to the numerical method for the p-biharmonic problem (and 
second order variational problems in general) is an appropriate definition of the 
Hessian of a piecewise smooth function. To formulate the general dG scheme for 
this problem from a variational perspective one must construct an appropriate 
notion of a Hessian of a piecewise smooth function. The finite element Hessian was 
first coined by |AM09) for use in the characterisation of discrete convex functions. 
Later in ILPllj was used in a method for nonvariational problems where the strong 
form of the PDE was approximated and finally put to use in the context of fully 
nonlinear problems in |LP12j . A generalisation of the finite element Hessian to 
incorporate the dG framework is given in |DP12j . which we also summarise here 
for completeness. 

Convergence of the method we propose is proved using the framework set out in 
[DPElOj where some extremely useful discrete functional analysis results are given. 
Here, the authors use the framework to prove convergence for a dG approximation 
to the steady state incompressible Navier-Stokes equations. A related but inde- 
pendent work containing similar results is given in (BO 09] where the authors study 
dG approximations to generic first order variational minimisation problems. 

The rest of the paper is set out as follows: The rest of this section introduces nec- 
essary notation and the model problem we consider. In ^we give some properties 
of the continuous p-biharmonic problem. In ^|3] we give the methodology for dis- 
cretisation of the model problem. In fj4]we detail solvability and the convergence of 
the discrete problem. Finally, in fjsjwe study the discrete problem computationally 
and summarise numerical experiments. 

Let i7 C R'* be a bounded domain with boundary dil. We begin by introducing 
the Sobolev spaces ^Cia78. .EvaQSj 



: J \(f>f < oo^ for p G [1, oo) and Loo(^^) = {0 : ess sup^ |0| < oo} , 

(1.1) 

W^in) = Lp{n) : D"(l) e Lp{n), for \a\ < 1} and := W^(l]), (1.2) 

which are equipped with the following norms and semi-norms: 

Mlm I IH' (1-3) 



l«llf.:=llx^lir,.,..- E llD"Hl£,(a) (1-4) 

|a|<fc 

\v\';^:^\vv:...,^.= E IP^^IImo) (1-5) 



|a|<fc 
|a|=A; 



ll^^ll? - M%{n) = ll«llw^(o)' (1-6) 

where a = {ai, . . . , ad} is a multi-index, \a.\ — ^f^^ ai and derivatives D° are 
understood in a weak sense. We pay particular attention to the cases I = 1,2 and 

Wl{n) := {0 e Wlin) : = (V0)^n = 0} . (1.7) 

In this paper we use the convention that the derivative Du of a function u : Q —> 
R is a row vector, while the gradient of u, Vu is the derivatives transpose, i.e., 
Vu = {Duy . We will make use of the slight abuse of notation, following a common 
practice, whereby the Hessian of u is denoted as D^m (instead of the correct VDu) 
and is represented by a d x d matrix. 
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Let L ~ L(a;, w, Vu, D^u) be the Lagrangian. We will let 



(1.8) 



be known as the action Junctional. For the p-biharmonic problem the action func- 
tional is given explicitly as 

1 



L(x. u. Vu, T> u) 



n 



\Au\' 



nP 



fu 



(1.9) 



where Aw :— trace(D^w) is the Laplacian and / G Lg(il) is a known source function. 

o o 

We then look to find a minimiser over the space Wp(r2), that is, to find u G Wp{il) 
such that 

Jf[u;p] = min ^[v;p]. (1.10) 

If we assume temporarily that we have access to a smooth minimiser, i.e., u € 
C^in), then, given that the Lagrangian is of second order, we have that the Euler- 
Lagrange equations are (in general) fourth order. 

Let X:Y = trace(X^l^) be the Frobenious inner product between matrices. We 
then let 



X 



then use 



dL 
d{X) 





-x\ .. 


■ -^1 








^d 
■ ■"'d. 




dL/dx\ 


.. dL/dxf 


_dL/dx\ 


. . dL/dx'^_ 



(1.11) 



The Euler-Lagrange equations for this problem then take the following form: 



D2 



dL 



dL 

du 



0. 



These can then be calculated to be 

^[u;p\ := A(\Auf-'^Au) = /. 



(1.12) 
i: 

(1.13) 
(1.14) 

Note that, for p — 2, the problem coincides with the biharmonic problem A^u = f 
which is well studied in the context of dG methods |SM07| iGNPOSl IGH09} e.g.]. 

2. Properties of the continuous problem 

To the authors knowledge the numerical method presented here is the first finite 
element method presented for the p-biharmonic problem. As such, we will state 
some simple properties of the problem which are well known for the problem's 
second order counterpart the p -Laplacian |Cia78[ fBEOSj . 

o 

2.1. Proposition (coercivity of Let u £ Wp(il) and f G Lg(il), where - + - = 

o 

1, we have that the action functional ]p\ is coercive over Wp(ri), that is, 

/[u-p]>C\u\lp-^, (2.1) 
for some C > and 7 > 0. Equivalently, let 



lAuF^^AuAv 



(2.2) 
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then we have that there exists a constant C > such that 

^{v,v;p)>C\v\l^ VveWli^}). (2.3) 

o 

Proof By definition of the Wp(il) norm we have that 

/[^■,P]^l\<2,p~f^- (2-4) 



Upon applying a Holder inequality we see 



1 

P 



(2.5) 



>^Kp-cil/llL,(n) 



The statement (2.3 1 is clear since 

^{v,v;p) = \v\l^, (2.6) 
thus concluding the proof. □ 

2.2. Proposition (convexity of L). The Lagrangian of the p~biharmonic problem 
is convex with respect to its fourth argument. 

Proof Using similar arguments to |Cia78l §5.3] (also found in [BL94]) the convexity 
of the functional J is a consequence of the convexity of the mapping 

^■.^eR^-m\P. (2.7) 
P 

□ 

2.3. Corollary (weak lower semicontinuity). The action functional ^ is weakly 

o 

lower semicontinuous over Wp{fl). That is, given a sequence of functions {%}je[M 

o 

who has a weak limit u G Wp(ri), then 

/[u-p\<\\mmi J[u,-p\. (2.8) 

Proof The proof of this is a straightforward extension of |Eva98[ §8.2 Thm 1] to 
second order Lagrangians, noting that ^ is coercive (from Proposition 2.11 and 
that L is convex with respect to its fourth variable (from Proposition 2.2 ). We omit 
the full details for brevity. □ 

2.4. Corollary (existence and uniquness). There exists a unique minimiser to the 
p-biharmonic equation. Equivalently there is a unique (weak) solution to the (weak) 

o 

Euler-Lagrange equations, find u G Wp{fl) such that 

[ lAuf^^AuAcj}^ [ f(j) \f<j)eWl{n). (2.9) 
Jn Jn 

Proof Again, the result can be deduced by extending the arguments in |Eva98| 



58.2] or (Cia78[ Thm 5.3.1], again, noting the results of Propositions 2.1 and 2.2 



the full argument is omitted for brevity. □ 

3. Discretisation 

Let ^ be a conforming, shape regular triangulation of 51, namely, is a, finite 
family of sets such that 

(1) K E ^ implies K is an open simplex (segment for d — 1, triangle for d — 2, 
tetrahedron for = 3), 
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(2) for any K,J & £^ we have that -RTfl J is a full subsimplex (i.e., it is either 
0, a vertex, an edge, a face, or the whole of K and J) of both K and J and 

(3) {jK^^K = n. 

The shape regularity of 5^ is defined as the number 

M(^) := inf ^, (3.1) 

where pk is the radius of the largest ball contained inside K and hx is the diameter 
of K. An indexed family of triangulations {.5^"}„ is called shape regular if 

/i:=inf/i(^") >0. (3.2) 

n 

We use the convention where /i : O IR denotes the meshsize function of ^ , i.e., 

/i(a;) := max/iif , (3.3) 

K3X 

which we shall commonly refer to as h. 

We let S be the skeleton (set of common interfaces) of the triangulation =^ and 
say e G if e is on the interior of Vl and e G dU, if e lies on the boundary dVt. 

We let P'^(^) denote the space of piecewise polynomials of degree k over the 
triangulation =^,i.e., 

P'=(^) = {0 such that <;i|ifeP'=(ii')} (3.4) 
and introduce the finite element space 

V:=DG(^,fc) = P'=(=^) (3.5) 
to be the usual space of discontinuous piecewise polynomial functions. 

3.1. Definition (finite clement sequence). A finite element sequence {u^.V} is a 
sequence of discrete objects, indexed by the mesh parameter h, individually rep- 
resented on a particular finite element space, V, which itself has discretisation pa- 
rameter h, that is, we have that V = V(/i). 

3.2. Definition (broken Sobolev spaces, trace spaces). We introduce the broken 

Sobolev space 

W^(^) := {(j) : (j)\K e W^(/r). for each K & 3r} . (3.6) 

We also make use of functions defined in these broken spaces restricted to the 
skeleton of the triangulation. This requires an appropriate trace space 

r(<f):= n L2(ai^)c n Wp-^(ir) (3.7) 

for p > 2, / > 1. 

3.3. Definition (jumps, averages and tensor jumps). We may define average, jump 
and tensor jump operators over T{S) for arbitrary scalar functions v e T{S) and 
vectors v S T{(o)'^. 

'^{v\k^ +v\kJ over ^ (3.8) 
v\dn on dCl 



Kt^ki+vkJ over^ (3.9) 
v\dn on 90 
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[wlxiW^i +w|x2"Jf2 over S (3.10) 
Xivn) Ian on do. 



\{v\KiVnK^ + {v\K^ynK.^ over S (3.11) 
[(v^n) \qq_ on 9r2 



u I-)- 



(8) n^i-j + ® riK^ over ff (3.12) 
(ti (g) n) laa on 9il 



We will often use the following Proposition which wc state in full for clarity but 
whose proof is merely using the identities in Definition |3.3| 



3.4. Proposition (clementwise integration). For a generic vector valued function 
p and scalar valued function (p we have 



/ div(p) (/)da; = ( — / p''' V /,,(/> da; + / (j)p^nxds 



(3.13) 



In particular, if we have p £ T{S U dVl)'^ and 4> £ T{S' U Sfi), the following identity 
holds 

I ^P'^K ds - / bl {{ 0} + / I p} ds = / Ip4>1 As. 

^ JdK Js Jsudn Jgudn 

(3.14) 



An equivalent tensor formulation of {3.13)~{3.14) is 



(3.15) 



/ D/iP(/)da;= ( ^ / P'XiV/i0da;+ / (g) n/^ ds 

In particular the following identity holds 
Y I ^P®riK ds^ [ Ipj^ l^} ds+[ {{ p} ds^ [ ds 

-^^ or JdK JS JSudn JSudn 

(3.16) 



The discrete problem we then propose is to minimise an appropriate discrete 
action functional, that is to seek u/i g V such that 

Jh[uh;p\ = inf /„K;p]. (3.17) 

3.5. Remark (motivation for discrete action functional). The choice of discrete 
action functional is crucial. A naive choice would be to take the piecewise gradient 
operators, substituting them directly into the Lagrangian, i.e., 

= / L{x,Uh,VhUh,^luh) . (3.18) 

This is, however, an inconsistent notion of the derivative operators (as noted in 
[BO09] i. 

Since for the biharmonic problem the Lagrangian is only dependant on the Hes- 
sian of the sought function, we need only construct an appropriate consistent notion 
of discrete Hessian. 

^^^^^^^ o 

3.6. Theorem (dG Hessian [DP12| ). Let v £ Wp(^) and suppose v and p are 
consistent numerical Ruxes representing approximations to v and Vu respectively 
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over the skeleton of the triangulation. Then the generaUsed dG Hessian, H[v] G 
V''^'' satisGes 



/ H[v] * = - / Vftw V,,$ - / 
Jn Jn Js 



sudn 
V<I> e V. 



m®{{pi+ j m 

sudn 



(3.19) 



Proof Note that, in view of Green's Theorem, for smooth functions, w e C^(il) n 
C^(il), we have 



Vu;«)V(/)+ / Vw®n</> y e C {fl) n C (fl) . (3.20) 

an 



As such for a broken function v G Wp(^) we introduce an auxiUary variable 
p = V/jW and consider the foUowing primal form for the representation of the 
Hessian of said function: For each K £ ^ 



/ H[v] $ = - / Vh$ + / n $ e V 

Jk Jk JdK 



p® q ^ — I V Dq 



K 



K 



q(g>nv y q eV , 



dK 



(3.21) 
(3.22) 



where =(D;i)^ is the elementwise spatial gradient. 
Noting the identity ( |3.16 1 and taking 

/ H[v] ^=Y. I * = E 

= - / Vh$+ / 

Jn JSudn 

Using the same argument for ( 3.22[ ) 

/ p®q^ E / P®^= E 

Jn Ke-Sr^ Kesr 

= - f vBhq + [ m<S)iq} + 





f / p (g) n 


^ Jk 


JdK 






Js 













(3.23) 



(3.24) 



m 



Note that, again making use of (3.161 we have for each q G and w S 

Hi(=^) that 



sudn 



{{'7}®M+ / I'ZLI^^}- (3.25) 



q <S) VhW = - / D/iQw; 

Taking w = w in (3.251 and substituting into (3.22) we see 

p®q^ I q®VhV+ I lv-vl®iq^+ I iv-vllql^^. (3.26) 
! Jn Jsvjdn Js 

Now choosing q — V?i<f> and substituting ( 3.26| ) into (3.21) concludes the proof. 

□ 

3.7. Example f |DP12p . An example of the possible choices of fluxes are 

« lovers (3.27) 

(3.28) 



[0 on dVL 
P=| Vftw} onS\Jdn. 
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The result is an interior penalty (IP) method |DD76| applied to represent the finite 
element Hessian 



Jn Jn JSudn 

(3.29) 



JSudfl 



3.8. Definition (lifting operators). From the IP-Hessian defined in Example 3.7 
we define the following lifting operator li,l2 : V ^ ydxd g^^j^ ^^j^g^^ 

/ /iK]$= / W®{{Vh$:& (3.30) 
Jn Jsudn 

I k[vh]'^ = - I IV,7.,1^|$}. (3.31) 
Jn Jsudn 

As such we may write the IP-Hessian as Jf : V ^ ydxri g^^.]^ ^j^g^^ 

/ H[vh]<^>= [ Dlvh<^> + li[vh]^ + hmvh + l2[vh]<i> V$eV. (3.32) 
Jn Jn 

4. Convergence 

4.1. Definition (dG norm). We define the dG norm for the p-biharmonic problem 
as 

\\vHr,a., w^'^-^^wim + h'^' wi^h-niwi^s) + h'-'' II wrL,(^) ^ (4-1) 

where A^- := trace(D^-) denotes the piecewise Laplacian operator. 

To prove convergence for the p-biharmonic equation we modify the arguments 
given in |DPE10j to our problem. To keep the exposition clear we will, where 
possible, use the same notation as in |DPE10] . 

4.2. Broken Sobolev embeddings. There is a wealth of material on embeddings 
for broken Sobolev spaces. We refer the reader to [LS031 iBOOOl IDPElOi e.g.]. We 
state some basic propositions, that is, a trace inequality and inverse inequality in 
Lp(ri), the proof of these is readily available in |Cia78[ e.g.]. Henceforth in this 
section and throughout the rest of the paper we will use C to denote an arbitrary 
positive constant which may depend upon and Q but is independent of h. 

4.3. Proposition (trace inequality). Let w/i G V be a finite element function then 
for p E (1, 00) there exists a constant C > such that 

\\vhk^iS')<Ch-'/P\\vHh^^n) (4-2) 

4.4. Proposition (inverse inequality). Let Vh G V be a finite element function then 
for p E (1, oo) there exists a constant C > such that 

\\ykV,\\lin)<Ch-n\vH\\lin) (4-3) 

4.5. Lemma (relating IHI^jq^ and IMI^^gt norms). For two integers s,t such that 
l<s<t<oowe have that there exists a constant C > such that 

\\vh\Us<C\\vHhG,t (4-4) 
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Proof The proof follows a similar line to [DPElOl Lem 6.1]. By definition of the 
norm we have that 



\dG,s 



IhMdG.s^ / \^kVHr + h'-' / \l^H.VHW + h'-'' / IKir- (4.5) 

Jn Jg Jg 

Now let us denote r ~ ^ and q = ^rrj, that is, we have that ^ + ^ = 1. Hence we 
may deduce that 



Jn Jg Jg 

< C'||v/J|rfG,t 

where we have used a Holder inequality together with 

1 - .s = 1 - ^ = i + and (4.7) 
1 - 2s = 1 - ^ = i + (4.8) 

and the shape regularity of ^ , concluding the proof. □ 

4.6. Definition (bounded variation). Let denote the variation functional de- 
fined as 



nu] ■■= sup { / wdiv0 : e [Clm^ ||</.||L^(a) < 1 



(4.9) 



The space of hounded variations denoted BV is the space of functions with bounded 
variation functional, 

BV := {(t> e Li(n) : r[</>] < c5o} . (4.10) 
Note that the variation functional defines a norm over BV , we set 

Mbv = nn]- (4.11) 

4.7. Proposition (control of the L^(17) norm |EGH10| ). Let u G BV then we 
have that there exists a constant C such that 

\\uh^^n)<C\\u\\sv (4.12) 

d-1 

4.8. Proposition (broken Poincarc inequality |BQ09| ). For f G V we have that 

||«^|lL,(i2) <c(£|V;,z;,,|+^|Kl|). (4.13) 

4.9. Lemma (control on the BV norm). We have that for each w/j e V and p £ 
[1, oo) that there exists a constant C > such that 

\\vh\\Bv<C\\vhhG,p (4.14) 
Proof Owing to [DPElOl Lem 6.2] we have that 

\\vh\\bv< I |V„«,.|+ / |M|. (4.15) 
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Applying the broken Poincare inequality given in Proposition 4.8 to the first term 
on the (4.15 ) gives 

\\vh\\sv<c([ |A„«,|+ / \lVHVkl\+ [ IWl) 



< c 

<C\\vh\\dG,j . 

Applying Lemma |4.5| concludes the proof. 



(4.16) 



□ 



4.10. Lemma (discrete Sobolev embeddings). For u/j G V there exists a constant 
C > such that 



(4.17) 



Proof The proof mimics that of the Gargliardo-Nirenberg-Sobolev inequality in 
|Eva98[ Thm 1, p.263]. 

We begin by noting that Proposition |4.7| together with Lemma |4.9| infers the 
result for p = 1, i.e., 

IkhllLuo) <C|I«/.Lg,i- (4-18) 
Now, we divide the remaining cases into two possibilities, p G (1, d) and p G [d, oo). 

Step 1. We begin with p G {l,d} and apply (4.181 to Vh = Iw/ip, where 7 > 1 
is to be chosen, obtaining 



\Wh 



< c 



/ liv 



h \Wh\ 



\Wh\ 



We proceed to bound each of these terms individually. 
Firstly, note that, by the chain rule 



7 \Wh 



\^h\wh\''\ = 
Using a Holder inequality it follows that 



where q = 



7-1 



(4.19) 



(4.20) 



(4.21) 



Now we must proceed to bound the skeletal terms appearing in (4.191. The jump 
terms here also act like derivatives in that they satisfy a 'chain rule' inequality, using 
the definition of the jump and average operators it holds that 



llV,Kn|< / 27SKr'SlV,u;,J 



< 27 



17-1 - 



(4.22) 



by a Holder inequality. 

Focusing our attention to the average term it holds that, by definition 



< 



Wh 



17-1 



h,{dK} 



(4.23) 
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11 



\Wh 



7-1 



\Wh 



|<j(7-l) 



(4.24) 



Choosing a = ^ such that the exponent of h vanishes and substituting into (4.22) 
gives 



(4.25) 



The final term is dealt with in much the same way. Again, using the 'chain rule' 
type inequality we see that 



/^-MiKrii<27 / h-'i\w,r'}\iwH 

<27 h^i\wur^ 



(4.26) 



which in view of (4.24) gives 



(/j-l""")' 









(4.27) 



again where a — - . 

Collecting the three bounds (4.21 ), (4.25) and (4.271 and substituting into (4.19) 
shows 



d-i 



\Wh 



|<j(7-l) 



\^hWh\ 



Lp(0) 



(4.28) 



The main idea of the proof is to now choose 7 such that = (7(7 — 1). Hence 
7 = . Using this and dividing through by the first term on the right hand 

side of (|428| yields 



d-l 1 



Now noting that 



yields 



d-l I _d-p 

d q dp 
h-^ = h^-P and 

Ik'Jlv(n) < W^hWdCp 



h 1 ^ hvh] 



(4.29) 



(4.30) 

(4.31) 
(4.32) 

(4.33) 



where p* — is the Sobolev conjugate of p. This yields the desired result since 
p* > p for p £ (1, d) and hence we may use the embedding Lp. (fi) CC Lp(f2). 
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Step 2. For the case p G [rf, oo) we set r — We note that r < d and that 

the Sobolev conjugate of r, r* = > r. Following the arguments given in Step 
1 we arrive at 

< \W\\iG,r- (4-34) 

Note that 



d-r d^-lf- 

d+p 



(4.35) 



Hence we see that 

W^hh^^n) = ll^''llL..(n) <C\\wh\\aG,r <C:\\wh\\aG,p' (4-36) 
where the final bound follows from Lemma |4.5[ concluding the proof. □ 



4.11. Theorem (stability). Let H[-] be defined as in Example 3.7 then the dG 
Hessian is stable in the sense that 



l^lvh - i?K]||Lp(0)<ixd < C[\\h[vh] + l2[vh]\\l^(^a)'ixa 



, . > (4-37) 



Consequently we have 

\\H[v,]\\l^,,^,^,<C\\v,r,G,p (4.38) 

Proof We begin by bounding each of the lifting operators individually. Let q = 
then by the definition of the Lp(ri) norm we have that 

IKiK]IIl,(^.)= sup / (4.39) 

Let Pv : L2(il) — > V denote the orthogonal projection operator then using the 
definition of (3.30) we see 



= sup / 



li[vh]Pvz 



)|V,.(Pv^)l 



= sup 

WllL,(^)ll«/i"V,(Pvz)}}||L^(^) 

< a sup r—r 

zeL,(f2) ll^llLg(O) 



< d sup 



(II/*-" WI1l,(<?)) (Ill h-^,,iFyz)}\\l 



(4.40) 



using a Holder inequality, followed by a discrete Holder inequality and where a € IR 
is some parameter to be chosen. 

Using the definition of the average operator we see 

||«/."v,(Pvz)}||^^(^)<i W'^^'^'^iPy^Km)- (4.41) 

Now using the trace inequality given in Proposition |4.3| we have 

III /i"V,(Pv ^) }||[^(^) <CY^ /i^-i II V,(Pv z)||[^(^) . (4.42) 
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Making use of the inverse inequality given in Proposition 4.4 
|||/i"V;,(Pvz 



9 



we see 

11"^ 



(4.43) 



We choose a = 2 — ^ such that the exponent of h in the final term of ( 4.43 ) is zero. 
Substituting this bound into (4.43) and making use of the stability of the L2(ri) 
orthogonal projection in hp{Q) we see that 



\\h[vH]\\l^n)<C 



Lp(<?) 



(4.44) 



< ch^-^p II 



The bound on l2[-] is achieved using much the same argument. Following the 
steps given in ( 4.40 ) it can be verified that 



Mvh]hrn)<:d^ sup 
zeL,(si) 



1/p 



1/9 



(4.45) 



L,(n) 



for some f3 G R. To bound the average term, we follow the same steps (without the 
inverse inequality) 



1 



< C h^^-' \\Py Z\\l 



(4.46) 



We choose /3 = 1 — ^ such that the exponent of h vanishes and substitute into 
^4A^ to find 



\h[vH]\\l^^^<C 



(4.47) 



<Ch'-P\\lv,Ml 



The result (4.37) follows noting the definition of H given in (3.32), a Minkowski 
inequality and the two results (4.44) and (4.47). 



To see (4.38) it suffices to again use a Minkowski inequality, together with (3.32) 



and the two results (4.44) and (4.47). 



□ 



4.12. Numerical minimisation problem and discrete Euler Lagrange equa- 
tions. The properties of the IP-Hessian allow us to define the following numerical 
scheme: To seek li/i e V such that 

[uh ; p] = inf [vh ; P] ■ (4.48) 
Let S^[vfi] '■— trace then the discrete action functional is given by 

JnP P ^Jgudn ' 

(4.49) 

where ct > is a penalisation parameter. 
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Let 



+ a 



Sudfl 



(4.50) 



The associated (weak) discrete Euler-Lagrange equations to the problem are to 
seek {uh, H[uh]) e V x V''^'' such that 



where H is defined in Example |3.7| 



(4.51) 



4.13. Theorem (coercivity). Let f G Lg(51) and {uh,^} be the finite element 
sequence satisfying the discrete minimisation problem {4.48) then we have that 
there exists constants C > and 7 > such that 



Equivalently let £/h{', SP) be defined as in \4.50) then 

£/h{uh,Uh;p) > C||u/i.|ldG,p- 
Proof We have by definition of IMI^g ^ that 



Uh\ 



IP 



(4.52) 
(4.53) 

(4.54) 
(4.55) 



We see by a Minkowski inequality that 

< W^hUH - + II^K]ll^,(o) 

Hence, using the stability of the discrete Hessian given in Theorem |4.11| we have 
that 

\\u,\Zg,p < m^'Mlm +(1 + c){h'-'' IIIv.MrL,(^) + h'-'" IIKlll^^(^)) 

<{l + C) £/hiuh,Uh;p) , 

(4.56) 



showing (4.53). The result (4.52) follows using a similar argument. 



□ 



4.14. Lemma (relative compactness). Let {w^, V} be a finite element sequence that 



is bounded in the 



\dG,p 



norm. Then the sequence is relatively compact in Lp{Q). 



Proof The proof is an application of Kolmogorov's Compactness Theorem noting 
the result of Lemma |4.10| which infers boundedness of the finite element sequence 
mLp{n). □ 

4.15. Assumption (approximability of the finite element space). Henceforth we 
will assume the finite element space V is chosen such that the L2(il) orthogonal 
projection operator satisfies: 

(4.57) 
(4.58) 
(4.59) 



lim ||«-Pvi;||l^(o) =0 



1™ II Vw-V,,(Pv 1^)111^^(12) 



and 



lim llw 



A choice of fc > 2 satisfies these assumptions. 
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4.16. Lemma (limit). Given a Gnite element sequence {vh,y} that is bounded in 

o 

the W'Wj^Q p norm, there exists a function v e Wp(il) such that as h we have, 
up to a subsequence, Vh ^ v weakly in hp{Q,). Moreover, H[vh] D^w weakly in 



Proof Lemma 4.14 infers that we may find a v € Lp(ri) which is the hmit of our 

o 

finite element sequence. To prove that v € Wp(i7) we must show that our sequence 
of discrete Hessians converge to D^u. 
Recall Theorem 4.11 gave us that 

\\H[v,M^^^^y,..<C\\vHha,p- (4.60) 

As such, we may infer the (matrix valued) finite element sequence {H[vh],'^'^^'^} 
is bounded in Lp(f2)''^''. Hence we have that H[vh] ^ X G Lp(r2)''^'' weakly for 
some matrix valued function X. 

Now we must verify that X — D^v. For each 4> G C[f'(ri) we have that 



Note that 



H[vn]<P= / Bivh<l>- / WhVhi^<l>+ / M®V</.. (4.61) 

(4.62) 



n 

n 



As such, we have that 



Xch^ lim / H\vi, 



lirn^ / VhT>^(l) (4.63) 



n 

'n 

and hence we have that X — D^w in the distributional sense. □ 

4.17. Lemma (apriori bound). Let f G Lg(il), with q = and let {u^,, V} be the 
Enite element sequence satisfying {4.48), then we have the following apriori bound: 

h/.LG,p<(C|l/llL,(o))'^'- (4.64) 

Proof Using the coercivity condition given in Theorem |4.13| and the definition of 
the weak Euler-Lagrange equations we have 

\\uh\\dGp ^ C£/h(uh,Uh;p) 

< C I fuh- 
Jn 

Now using a Holder inequality and the discrete Sobolev embedding given in Lemma 
14.101 we see 

iI"/.ii!;g,p <G|i/iiL,(o)ii"'^iiL,(a) 

<C\\fh^^n)\Wd,G,,- ■ 
Upon simplifying, we obtain the desired result. □ 

4.18. Theorem (convergence). Let / G Lq(17), with q = and suppose {uh,^} 
is the finite element sequence generated by solving the nonlinear system {4.51}, 
then we have that 
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• Uh ^ u in Lp(r2) and 

• H[uh] D^M in Lp{n)'^'"^. 

where u € Wp(ri) be the unique solution to the p-biharmonic problem (1.14). 

Proof Given / £ Lg(f2) we have that, in view of Lemma 4.17 the finite element 
sequence {uf^y} is bounded in the norm. As such we may apply Lemma 



4.16 which shows that there exists a (weak) hmit to the finite element sequence 
{uh, V} which we shall call u* . We must now show that u* = u, the solution of the 
p-biharmonic problem. 

^[■] is weakly lower semicontinuous, hence we have that 

"1 



By Corollary 



2.3 



< liminf 

< liminf 



\S![uh_ 



MP 

llLp(O) 



Lp(n) 



fUh 



(4.67) 



hUh_ 



IP 



hm inf Jfh [uh] 



Now owing to Assumption 4.15 we have that for any v e C[^(51) that 



lim inf 

/i-s-O 



'"-{h^~^mH{Vyv)l\\l^ 



h'-^p I 



pvt'irL,(o) 



(4.68) 



hm inf Pv v] 



By the definition of the discrete scheme we have that 

^K]< AK]< A[Pv«] = 



(4.69) 



Now, since v was a generic element we may use the density of C[5°(f^) in Wp(ri) 
and that since u is the unique minimiser we must have that u* ~ u. □ 

4.19. Remark (provable rates for the 2-biharmonic problem). In the papers |SM07[ 
IGH09j rates of convergence are given for the 2-biharmonic problem, these are 

\\u - UhW = 0(/i2) for /c = 2 ||w - UhW = 0{h^+^) for /o 2 (4.70) 

II"-"'.Lg,p = 0('^'"')- (4-71) 
Note that for piecewise quadratic finite elements the convergence rate is suboptimal 
in U{^)- 

5. Numerical experiments 

In this section we summarise some numerical experiments conducted in the 
method presented in ^J3] 

5.1. Remark (implementation issues). The numerical experiments were conducted 
using the DOLFIN interface for FEniCS |LW10| . The graphics were generated using 
Gnuplot and ParaView . 

For computational efficiency, we chose to represent 2![uh\ as an auxiliary variable 
in the mixed formulation, which only requires one additional variable, as apposed 
to the full discrete Hessian H[u}-^ which would require (P (or <^ if one uses 
symmetry of H). We note that this is only possible due to the st;ructure of the 
problem, i.e., that L = i(a;, u, Vw, Aw) and would not be possible in a general 
setting. 
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5.2. Benchmarking. The aims of tliis section are to test the robustness of the 
numerical method for a model test solution of the p-biharmonic problem. We show 
the method achieves the provable rates for p ^ 2 (Figure [l]) and numerically gauge 
the convergence rates for p > 2 (Figures [2] and [s]). To that end, we fix d = 2, let 
X ={x,yy and choose / such that 



u{x) :— sin (27ra;)^ sin (27r?/)^ . 



Note that this is comparable to the numerical experiment |GH091 §6.1]. 



(5.1) 



Figure 1. §5.2| - Numerical experiment benchmarking the nu- 
merical method for the 2-biharmonic problem. We fix / such that 



the solution u is given by (5.1 1. We plot the log of the error to- 
gether with its estimated order of convergence. We study the Lp{V,) 
norms of the error of the finite element solution Uh as well as the 
represented auxiliary variable ^[m/i] for the dG method ( 4.51| ) with 
k = 2,3,4. We also give a solution plot. We observe that the 
method achieves the rates given in Remark |4.19| 
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(a) finite element approximation to | |5.1[ l. 



(b) k = 2, piecewise quadratic FEs. 
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10000 100000 
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(c) fc = 3, piecewise cubic FEs. 



(d) = 4, piecewise quartic FEs. 



5.3. Remark (computational observations). During the experiments we noted a 
superconvergence for both the Lp{n) errors of the solution as well as the auxiliary 
variable for odd values of p and even values of k. Computationally we observed 
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Figure 2. §5.2| - The same test as in Figure [T] for the 3- 
biharmonic problem, i.e., p = 3. In addition, we test sextic el- 
ements. Note there is a superconvergence for quadratic, quartic 
and sextic elements, see Remark 15.31 
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(a) k = 2, piecewise quadratic EEs. 



(b) k = 3, piecewise cubic EEs. 
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(c) k = 4, piecewise quartic EEs. (d) k = 6, piecewise sextic EEs. Note that in 

this experiment numerical precision becomes an 
issue when studying the ljp(Q) error of the so- 
lution, but superconvergence can still be seen in 
the error of the auxiliary variable. 



that 



\U - Uh\ 



Lp(n) 



0{h P ) when fc > 2, fc is even and p is odd 
0(/i'^+^) otherwise 



when k 
1 



and that 



I At 



Lp(n) 



k-l + 



when k is even and p is odd 



(5.2) 



(5.3) 



0{h'' ^) otherwise . 

Compare with Figures |3(a)| [3(c)] and [3(d)] and Figures [5(a)] [5(c)] and [5(d)] 

5.4. Remark (representation of H). Note that the dG Hessian H may be repre- 
sented in a finite ele ment space with different degree to Uh S V. Let W := P'^^^{£^), 
the proof of Theorem 3.6 infers that we may choose to represent H[uh] G W''^''. For 
clarity of exposition we choose to use H[uh] E \dxd^ however, we see no difficulty 
extending the arguments presented to the lower degree dG Hessian. 
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Figure 3. §5.2| - The same test as in Figure [2] for the 4- 
biharmonic problem, i.e., p = A. 
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(a) k = 2, piecewisc quadratic FEs. 

^ ^ ^ IDelta u-tetiV^lIlt^ ~ 



(b) k = 3, piecewise cubic FEs. 



EOC = 2.99 

-EOC = 3.00 




^EOC = 5.00 
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(c) fc = 4, piecewise quartic FEs. 



(d) A: = 6, piecewise sextic FEs. 



Numerically we observe the same convergence rates as in Remark 5.3 for the 
lower degree dG Hessian. 

6. Conclusion and outlook 

In this work we presented a dG finite element method for the p-biharmonic 
problem. To do this we introduced a auxiliary variable, the finite element Hessian 
and constructed a discrete variational problem. 

We proved that the numerical solution of this discrete variational problem con- 
verges to the extrema of the continuous problem and that the finite element Hessian 
converges to the Hessian of the continuous extrema. 

We formulated the problem in this way due to the applications to discrete sym- 
metries, such as discrete versions of Noether's theorem, which will be studied in a 
forthcoming publication. 

We foresee that this framework will prove useful when studying other (possibly 
more complicated) second order variational problems, such as discrete curvature 
problems like the affine maximal surface equation, which is the topic of ongoing 
research. 
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